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We construct a class of projected entangled pair states (PEPS) which is exactly the resonating 
valence bond (RVB) wavefunctions endowed with both short range and long range valence bonds. 
With an energetically preferred RVB pattern, the wavefunction is simplified to live in a one param- 
eter variational space. We tune this variational parameter to minimize the energy for the frustrated 
spin 1/2 Ji — J2 antiferromagnetic Heisenberg model on the square lattice. Taking a cylindrical 
geometry, we are able to construct four topological sectors with even or odd number of fluxes pen- 
etrating the cylinder and even or odd number of spinons on the boundary. The energy splitting 
in different topological sectors is exponentially small with the cylinder perimeter. We find a power 
law decay of the dimer correlation function on a torus, and a InL correction to the entanglement 
entropy, indicating a gapless spin liquid phase at the optimum parameter. 

PACS numbers: 75.10.Kt,75.10. Jm 



Introduction - Resonant valence bond (RVB) states, 
which were first proposed by Anderson [JJ to describe a 
possible ground state for the 5 = 1/2 antiferromagnetic 
Heisenberg model on a triangular lattice, and later to ex- 
plain the possible mechanism of high-T c cuprates (2] [3] , 
provide us a rich tool box to construct the so called spin 
liquid states. Rokhsar Kivelson (RK) wavefunction [3], 
which is an equal weight superposition of nearest neigh- 
bor (NN) dimer coverings, is a critical spin liquid state [3] 
on square lattice; whereas a gaped Z2 spin liquid state 
on kagome and triangular lattice [6] [7] . The equal weight 
superposition of the NN RVB state on square lattice was 
shown to be critical [5J |S] • Several numerical work [TTJ1 - 
IT3"] have demonstrated that the equal weight NN RVB 
states on the kagome and triangular lattices are Z2 spin 
liquid states. 

Recently numerical breakthroughs claimed a spin liq- 
uid ground state for the Kagome Heisenberg model [HI 
H5] and the frustrated spin 1/2 J\ — J2 antiferromagnetic 
(AF) Heisenberg model on the square lattice [TH1 H7] - 
However, these work did not give direct access to the 
topological nature of the spin liquid states, therefore, a 
simple variational wavefunction approach is highly desir- 
able. Although the variational energy of the NN RVB 
state on the kagome lattice [JTJ [13] is still higher than 
the energy obtained via the density matrix renormal- 
ization group (DMRG) method [2], the topological na- 
ture is well understood within the formalism of the pro- 
jected entangled pair states (PEPSs) [IT] . On the other 
hand, from a projective wavefunction [18 approach sup- 
plemented by a projective symmetry group (PSG) anal- 



ysis all possible spin- liquid states on the triangular [15] . 
Kagome 19-22 and Honeycomb [53] lattices have been 
obtained and classified but, for all lattices, the energeti- 
cally favorable states are believed to involve longer range 
RVB. As a result, it is natural to think that a general 
RVB state within the PEPS formalism is a more practical 
variational wavefunction, where one can gain simultane- 
ously an optimized energy and a comprehensive picture 
of the topological nature. 

In this paper, we introduce a general RVB state writ- 
ten as a D = 3 PEPS, different from Ref. [HH EE], i.e. 
it includes valence bonds of all length (although with 
a bond amplitude decaying exponentially fast with the 
bond length). With a properly chosen singlet sign con- 
vention that meets all lattice symmetries on the square 
lattice, we minimize the energy of the spin 1/2 J\ — J 2 
AF Heisenberg model at J2 — 0.5 Ji against a single vari- 
ational parameter c governing the decay amplitude of the 
long range valence bonds. The idea is therefore to intro- 
duce a simple yet competing wavefunction that enables 
us to fully understand the topological properties of the 
ground state of the frustrated magnets. 

RVB states in PEPS formalism - The equal weight 
superposition of the NN RVB states can be constructed 
using a PEPS with bond dimension D = 3 as following: 
each physical site has 4 virtual spins attached, each of 
which spans a virtual dimension of spin 1/2 © 0. From 
the bond point of view, every pair of the NN virtual spins 
is projected to a block diagonal virtual spin singlet state: 

|5) = |Q1)-|10) + |22), (1) 



2 



(a) 



(a2) 



(b) I (c) ^ (d)l 

(bl) . (cl) (dl) . 



(e) 



(b2). (c2) (d2K 



FIG. 1: (a-d) Red thick lines denote mappings of a virtual 
spin 1/2 to a physical spin 1/2, the blue dash lines represent 
the singlet pairings of two virtual spin-l/2s, and the black 
thin lines represent the virtual spin states, (e) The local 
sign convention for the bond singlets (in green arrows) and 
the corner singlets (in blue arrows). 



here the virtual indices "0,1" span the subspace of spin 
1/2 and virtual index "2" spans the subspace of spin 
0. At the physical site, a projector enforces one of the 
virtual spins with its spin 1/2 subspace to be mapped to 
the physical spin 1/2 state and the rest of virtual spins 
to stay in the spin subspace, i.e. the "2" state, 



= YJ\ t)<0|* + HXMk) ® (222| /fc , 



(2) 
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here subscript "/fc" stands for all except k. This PEPS, 
by contracting the virtual index of each S at the bond 
and each V\ at the vertex, represents exactly the equal 
weight NN RVB states. 

To allow long distance singlet pairings, we need spins 
to teleportate: enforcing a singlet between site i and j 
that are already paired in singlets (s±,i) and {s%, j) will 
generate a singlet pair (si,S2). The following projector 
realizes spin teleportation without increasing the bond 
dimension, 
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(I t)<0|i + I 4>>(1| ( ) ® (2\j ® (e\ kl , (3) 



here |e)jw = |01)jm — |10)jy, and it forces spins connected 
via this site by bonds k and I into a singlet. A general 
RVB wavefunction is a parameter c weighted combination 
of projectors V = V1+CV2 at each vertex V = {v} traced 
out with the bond singlets S at each bond B = {&}, 
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The sign convention and symmetries - Fig. [ij a-d) 
enumerates 4 possible V\ projectors and 8 V2 projectors 
at each vertex. The bond singlet S is chosen such that 
NN singlets point from sublattice A to B; the corner sin- 
glets, which plays the role of singlet teleportation, are ori- 
ented counter clockwise and preserve all lattice symme- 
tries. The sign convention is demonstrated in Fig. [l|e). 

The NNN singlet arises through two bonds singlets and 
one corner singlet, as in Fig. [5Ja). However the weight of 
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FIG. 2: (a) The NNN singlet (si,s 2 ) pairs through 2 bond 
singlets and 1 corner singlet via 2 paths, which cancel each 
other, (b) The next allowed AB sublattice pairing (si,S2) 
through 3 bond singlets and 2 corner singlets, (c) The weight 
distribution in the Liang-Doucot- Anderson picture with a lin- 
ear color scale for a 8 x 8 torus at c = 0.35. h(dx, dy) is the 
weight of a singlet at separation (dx,dy), where fo(l,0) = 1. 
The plotted color scale takes the square root of h(dx,dy) to 
magnify the weight of the long range singlet. 



a diagonal singlet is comprised of two shortest paths of 
equal magnitude but opposite sign, thus the net weight 
of the diagonal singlet is zero. The only shortest path 
to build the next range AB sublattice singlet is shown in 
Fig. [2^b) , and it consists of three bond singlets and two 
corners. The sign of the next range AB singlet is pointing 
from sublattice A to B. In general, no A A pairings survive 
(see Appendix) and all AB pairings point from sublattice 
A to B. To verify this result, we implement a Monte 
Carlo (MC) sampling of the singlet distribution of Q 
and calculate the weight h(dx, dy) defined in Ref. [18] 
as a function of separation. The result is presented in 
Fig. [2jc) and is consistent with the above analysis. 

String picture - Since the RVB bonds only connect 
sites on different sublattices, we can view such a RVB 
state as a liquid state of oriented strings. Indeed, choos- 
ing a reference VB configuration, any VB configuration 
can be viewed as a closed oriented string configuration: 
the RVB bonds in the VB configuration are regarded as 
a piece of the closed string pointing from the A to the 
B sublattice, while the RVB bonds in the reference VBo 
configuration are regarded as the complementary piece 
pointing from the B to the A sublattice. 

If such a superposition of closed orientable string states 
indeed represent a liquid state of closed strings, then 
the entanglement entropy for a such a state in a region 
A has the form Sa = clLa — \\q.La + b, where La is 
the length of the perimeter of region A. To understand 
such a result, we view a string as a flux line and the 
close-string condition implies that the flux is conserved. 
Therefore, the total flux going through the perimeter of 
A is zero. If we had ignored the flux conservation and 
assumed that the flux could fluctuate freely and inde- 
pendently, then the entanglement entropy would have an 
exact area-law form Sa = o,La and the typical amount 
of flux through the perimeter would be proportional to 
y/T/A. So if we restrict the amount of flux through the 
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Z = diag(l,l,-l) 



FIG. 3: (a) The parity quantum number G v (Gh) is defined by 
counting the number of singiets (blue lines) modulo 2 crossed 
by a vertical loop (horizontal line) in dashed black line, (b) 
By cyclically permuting all spins in a loop winding around 
the cylinder, a Gh = 1 configuration in (a) is transformed 
into a Gh = — 1 configuration in (b). Singlets in green lines 
are affected by the permutation operator, (c) An odd-flux 
state is defined by inserting operator Z on all vertical bonds 
crossed by a horizontal line to Q. 



perimeter to be zero, the entanglement entropy will be 
Sa = clLa — \ii\/La + b = clLa — \\x\La + b. The IuLa 
dependence in the entanglement entropy implies that the 
liquid of orientable closed strings must be gapless. This 
property is confirmed by our numerical calculation. 

Variational ground state energies at J 2 = 0.5Ji 
- We consider the general RVB wavefunction on a cylin- 
der with finite cylindrical circumference N v and infinite 
horizontal length Nh = 00. The physical properties are 
determined by the eigenvector with the largest eigenvalue 
of the transfer matrix. Let us introduce a horizontal 
(vertical) parity number Gh (G v ) which is defined by 
counting the number of singlets modulo 2 that cross a 
horizontal (vertical) line joining the two boundaries of 
the cylinder (going around the cylinder). The two states 
with Gh = ±1 (+ is even and — is odd) are orthogonal to 
each other in the thermodynamic limit, and can be trans- 
formed from one to another by a cyclic spin permutation 
LTq operator winding around the cylinder as illustrated 
in Fig. [^b). 

The Gh = ±1 states are not the minimally entangled 
states (MESs) 24 . However their superpositions, 



|*(±)) = |*)< 



(5) 



with a relative ± sign are. A good reason for it is that 
these states ^ can be written as simple PEPSs: |^(+)) 
state is the state corresponding to Q, and |4 r (— )) state 
is obtained by inserting a "vison" line to the PEPS for 
state | V E'(+)) (by putting Z — diag(l, 1, — 1) operators 
on the vertical bonds crossed by a horizontal line), as 
in Fig. [3jc) . These MESs are referred to as the even- 
flux (!#(+))) states and the odd-flux (!*(-))) states 
which stands for even and odd number of flux penetrat- 
ing the cylinder. We will show next that these four states 
|\£(±l)) e / are (bulk) ground states of the gapless spin 
liquid state. 
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FIG. 4: Ground state energies (per site) computed for the 
Ji — J2 AF Heisenberg model at J2 = 0.5Ji, as a function of 
parameter c for the e/o topological sectors of (a) a cylinder 
(N v = 4,6,8) with even-flux (N v /2 even) or odd-flux (N v /2 
odd), and (b) a strip geometry with N v — 10, ■ • ■ , 18. For 
both cases, energy minimizes at c = 0.35(1). 



The variational energies of |\&(±)) e / on cylinders with 
finite perimeter N v = 4, 6, 8 (N v must be even, otherwise 
the system dimerizes) are computed exactly via the trans- 
fer matrix method and shown in Fig.[4|a) (an even or odd 
number of flux is chosen to provide the lowest energy) as 
a function of the variational parameter c. The best vari- 
ational energy for the spin 1/2 Ji — J 2 AF Heisenberg 
model at J 2 = 0.5 Ji is c = 0.35(1) with N v = 4,6,8. 
To access larger system size, we study a complemen- 
tary geometry where the cylinders are cut open, with 
the top and bottom vertical virtual spins set to "2"s. 
We call them the finite (N v ) width strips. For a con- 
tractible geometry as strips, the flux parity is no longer 
meaningful, but the boundary parity quantum number 
G v — e(o) still holds. We simulate the leading eigenvec- 
tors of the transfer matrix of the strips by the matrix 
product states (MPS) with the same quantum number 
G v in both bra and ket. The ground state energies as a 
function of the variational parameter c for both sectors 
(| 4") e/o) of the strips with N v = 10, ■ • • , 18 are presented 
in Fig. |4|b). The best variational energy for J 2 = 0.5Ji 
is at c = 0.35, which is in good consistency with the case 
of finite perimeter cylinders. 

The variational energies of the even and odd sectors at 
the optimum parameter c = 0.35 as a function of inverse 
width 1/N V are shown in Fig. [5ja,b), with cylinders of 
size N v = 4, 6, 8 and strips up to N v = 36. A linear re- 
gression is applied to the even sector of the strips and a 
thermodynamic limit of — —0.48612(2) is obtained. 
This energy is competing on the third decimal digit to 
the best variational estimate of = —0.4943(7) with 
a D = 9 PEPS [17] , let alone the fact that here we vary 
only one variational parameter in a D = 3 PEPS. A con- 
jecture about the ground state energies of the gapless and 
gaped spin liquid states is that the energy splittings be- 
tween different topological sectors become exponentially 
small with the system size. This conjecture is verified 
in Fig. [5jc,d) presenting on semi-log scales the energy 
difference between all sectors and for cylinders and 
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FIG. 7: Renyi entanglement entropy &(£) for an area L x L 
on a torus of size 2L x L with L = 4, • • • , 20. The fitted red 
line is of form S 2 (L) = a x L - flnL + 6i for L € [6 : 18]; 
the fitted green line is 5*2 (L) = a 2 L + b 2 for L £ [6 : 14] and 
b 2 = -0.68(1). 



FIG. 5: Ground state energies at c = 0.35 as a function of 
1/N V for all topological sectors in (a) cylinders and (b) strips. 
The extrapolated ground state energy from the even sector 
of strips is — —0.48612(2). The ground state energy 
splitting between the lowest energy sector and other sectors 
as a function of N v for (c) cylinders and (d) strips. The 
splitting vanishes exponentially with size N v . 
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FIG. 6: The dimer (a) and the spin (b) correlation function 
for all topological sectors on a jV„ = 8 cylinder at c = 0.35. 
The dimer (c) and the spin (d) correlation function on a torus 
with L — 8, 16, 32, 64 at c = 0.35. Note that different scales, 
log-log in (c) and semi-log otherwise, are used. 



between the two existing energy sectors for strips. 

Correlation functions and entanglement en- 
tropy - We define the spin and dimer correlation func- 
tions as the ground state expectation values C(r) — (S • 
S r )and£>*(r) = ((So-S 1 )(S r -S r+1 )-(S -S 1 )(S r _ r S r )>. 
Fig. [6^ a) plots the dimer correlation functions on a cylin- 
der with N v — 8 for all topological sectors at the optimal 
parameter c = 0.35. The odd sectors have very slowly 
decaying dimer correlations due to an odd number of 
spinons sitting on the boundaries, thus the system effec- 
tively becomes an odd-width cylinder and the Majumdar- 
Ghosh kind of degeneracy emerges. We can eliminate 



the boundary effect by setting the system on a torus and 
carrying a variational MC simulation for PEPSs [3S] . We 
found the dimer correlation function exhibits a power law 
decay D*(r) ~ with a — 1.4, as shown in Fig. [^c). 

In contrast, the decay of the spin correlation function for 
all sectors on a cylinder or on a torus remains exponen- 
tial with a correlation length £ s 1.1, as evidenced in 
Fig. [6]jb,d). Fig. [T] shows the Renyi entropy 5*2 (L) of an 
area L x L on a 2L x L torus for size L — 4, 6, • • • , 20. 
The fitted line of aL— ^\nL + b reflects the InL correction 
from the oriented string picture. The simulation is done 
via MC sampling of the RVB configuration . Finally, 
we also would like to point out that the logarithmic cor- 
rection is very hard to be detected on small system size. 
If we fit S2 with a form aL + b on small system size, we 
find that the constant b — —0.68(1), which is very close 
to — ln2. Such an observation implies that the observed 
— In2 constant in DMRG calculation [TB] is insufficient to 
rule out the possibility of gapless spin liquid ground state 
for Ji — J2 Hcisenberg model on square lattice. 

Conclusion and outlook - We constructed a class of 
projected entangled pair states which exactly represent 
general RVB wavefunctions with all bond length contri- 
butions. Upon choosing an energetically preferred RVB 
pattern, we are able to build a one-parameter manifold 
of variational RVB D = 3 PEPSs which preserve all lat- 
tice symmetries. Minimization of the variational energy 
for the frustrated spin 1/2 J\ — J 2 Antifcrromagnetic 
(AF) Heisenberg model on the square lattice yields, at 
J2 = 0.5Ji, an energy Erya = —0.48612(2) per site in the 
thermodynamic limit. In the case of a cylinder geom- 
etry, four orthogonal topological states were identified, 
namely the even-flux and odd-flux states with even and 
odd number of spinons on the boundary. We found the 
dimer correlation function decays algebraically while the 
spin correlation function still decays exponentially. The 
entanglement entropy scaling reveals InL correction to 
the area law. Both features point towards the gapless 
spin liquid nature of our constructed RVB wavefunction. 

Previous valence bond MC simulations have proposed 
wavefunctions which violate the Marshall's sign rule by 
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a single negative pairing magnitude h(2, 1) [571 HE], how- 
ever our PEPS wavefunction constructed to meet a neg- 
ative /i(2, 1) condition does not gain an optimized energy 
except on a very small 4x4 torus. 

The PEPSs construction of the general RVB states 
can be applied to other bipartite and non-bipartite lat- 
tices, where the Schwinger boson spin liquid states under 
the projective symmetry group (PSG) analysis have been 
found [EH231 121] , but for which thermodynamic energies 
and correlation functions are still unknown due to a neg- 
ative sign problem in the valence bond MC simulations. 
Within the PEPS formalism all of these can be easily 
studied. Our PEPS construction of the RVB states can 
be further generalized to accommodate more complicated 
pairing pattern which can improve further the ground 
state energy although possibly requiring a larger bond 
dimension. 
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FIG. 8: (Color online) For any teleportation path contributing to the AA (BB) pairing, there always exists a dual teleportation 
path with opposite sign. Therefore the AA {BB) pairing amplitudes are exactly zero in |5 , }rvb- 
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Supplementary material 

'J / )rvb has no AA (BB) pairing 

In the main text, we have illustrated why the leading order contribution of NNN pairing vanishes in I^rvb- Here 
we would like to show that, indeed, all the AA (BB) pairing vanish in |^)rvb exactly up to any order in c. As 
explained in the text, the pairing amplitude between two sites i and j is given by the sum of all teleportation paths 
that connect i and j. As seen in Fig. |5J the solid red line a is a typical teleportation path that contributes to the 
pairing (ij). According to the local rule of the teleportation path (each path must turn at a vertex), it is not hard to 
see that any teleportation path contributing to the AA (BB) pairing must go through an odd number of corners and 
an even number of links. On the other hand, for a given teleportation path a, we can always find a dual teleportation 
path b which goes though the same number of corners and links. For any pair of dual a and b paths, the sign 
contribution from all the links (corners) are the same (opposite). Therefore, the total contribution from the pair of 
teleportation paths a and b always vanishes. Thus, we have proved that |\P)rvb h as n0 (BB) pairing. 



Finite size extrapolation of the optimum parameter c 



We provide here a more accurate determination of the optimum parameter c. For this purpose, we use a quadratic 
function to fit the even sector (|4 r ) e ) energy of a finite width (N v ) strip as a function of parameter c (around 
the minimum), and we extract the optimum parameter c op t(N v ) and the minimum energy iJ m i n (A^„) for N v = 
10, • • • , 24. The results are presented in Fig. ^a). Taking the optimum parameter c opt (N v ) and extrapolating it to 
the thermodynamic limit as a function of inverse width 1/N V , we find c opt (oo) = 0.356(1) (see Fig. [9](b)). Again, 
taking the minimum finite width energy E m - la (N v ) and extrapolating it to the thermodynamic limit as a function of 
l/N v , we obtain a thermodynamic energy = —0.48620(1 ) (see Fig.[9^c)). 
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FIG. 9: (Color online)(a) A quadratic fit of the even sector energy E(N V ) of a strip as a function of c around its minimum with 
N v — 10, • • • , 24. (b) Optimum parameter c op t(N v ) plotted as a function of 1/N V . A linear regression gives c opt (oo) = 0.356(1). 
(c) Minimum energy i? m i n (A„) plotted as a function of 1/N V . A linear regression gives E(oo) = —0.48620(1). 



